package fem2.material.cohesive_law;

import fem2.Element;
import fem2.NotImplementedError;
import fem2.Variable;
import fem2.element.crack_sda.ConstantJumpSDAElement;
import fem2.element.crack_sda.SDAElement;
import fem2.enu.VariableName;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import math2.MutableDouble;

/**
 * Mode I cohesive law based on Sancho paper: An embedded crack model for finite
 * element analysis of concrete fracture.<br>
 * Note: this class works with constant opening only
 * 
 * @author hbui
 * 
 */
public abstract class Mode1SanchoCohesiveLawImplex extends CohesiveLaw {

	private double H;

	public Mode1SanchoCohesiveLawImplex(double Gf, double ft) {
		super(Gf, ft);
		H = -Math.pow(ft, 2) / (2 * Gf);
	}

	/**
	 * compute value of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected abstract double computeSofteningLaw(double x);

	/**
	 * compute value of derivative of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected abstract double computeSofteningLawDerivative(double x);

	@Override
	public IVector computeJumpTraction(Element e, IVector w) {
		IVector t;

		// /*
		// * compute traction by checking for loading / unloading
		// */
		// double wtilde_alg = getVariable(e).getScalarVariable(
		// VariableName.EQUIVALENT_JUMP_ALGORITHMIC).getValue();
		//
		// /*
		// * check for onset of opening, using linear cohesive law instead
		// */
		// if (wtilde_alg < MathUtilities.tol6) {
		// double ft = computeSofteningLaw(0);
		// int dim = e.getMeshPart().getDimension();
		// t = new ArrayVector(dim);
		// t.set(0, ft);
		// return t;
		// }
		//
		// t = (IVector) w.clone();
		// double wtilde_old =
		// getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP, -1);
		// double wtilde_old_old =
		// getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP, -2);
		// double delta_wtilde = wtilde_old - wtilde_old_old;
		// double f_old =
		// getVariable(e).getScalarVariable(VariableName.STRESS_LIKE_VARIABLE,
		// -1);
		// double H = computeSofteningLawDerivative(wtilde_old);
		// double f = f_old + H * delta_wtilde;
		// BLAM.scale(f / wtilde_alg, t);

		int dim = e.getMeshPart().getDimension();
		t = new ArrayVector(dim);
		t.set(0, ft + H * w.get(0));

		return t;
	}

	@Override
	public IMatrix computeJumpStiffness(Element e, IVector w) {
		int dim = e.getMeshPart().getDimension();
		IMatrix C;

		// double wtilde_alg = getVariable(e).getScalarVariable(
		// VariableName.EQUIVALENT_JUMP_ALGORITHMIC).getValue();
		//
		// /*
		// * check for onset of opening, using linear cohesive law instead
		// */
		// if (wtilde_alg < MathUtilities.tol6) {
		// C = new Array2DMatrix(dim, dim);
		// C.set(0, 0, this.H);
		// return C;
		// }
		//
		// C = new IdentityMatrix(dim);
		// double wtilde_old =
		// getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP, -1);
		// double wtilde_old_old =
		// getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP, -2);
		// double delta_wtilde = wtilde_old - wtilde_old_old;
		// double f_old =
		// getVariable(e).getScalarVariable(VariableName.STRESS_LIKE_VARIABLE,
		// -1);
		// double H = computeSofteningLawDerivative(wtilde_old);
		// double f = f_old + H * delta_wtilde;
		// BLAM.scale(f / wtilde_alg, C);

		C = new Array2DMatrix(dim, dim);
		C.set(0, 0, this.H);

		return C;
	}

	@Override
	public void Initialize(Element e) {
		Variable v = new Variable(2);
		v.registerScalarVariable(VariableName.EQUIVALENT_JUMP, true);
		v.registerScalarVariable(VariableName.STRESS_LIKE_VARIABLE, true);
		v.registerScalarVariable(VariableName.EQUIVALENT_JUMP_ALGORITHMIC, false);
		setVariable(e, v);
	}

	@Override
	public void InitializeSolutionStep(Element e) {
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {

		/*
		 * explicitly extrapolate the equivalent jump
		 */
		double wtilde_old = getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP, -1);
		double wtilde_old_old = getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP, -2);
		double delta_wtilde = wtilde_old - wtilde_old_old;
		double wtilde_alg = wtilde_old + delta_wtilde;
		getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP_ALGORITHMIC).setValue(
				wtilde_alg);

		double wtilde;
		if (e instanceof ConstantJumpSDAElement) {
			IVector alpha = ((SDAElement) e).getCrackOpeningParameter();
			wtilde = BLAM.norm(BLAM.NORM_TWO, alpha);
		} else {
			throw new NotImplementedError();
		}

		/*
		 * update equivalent jump implicitly
		 */
		MutableDouble wt = getVariable(e).getScalarVariable(VariableName.EQUIVALENT_JUMP);
		if (wt.getValue() < wtilde) {
			wt.setValue(wtilde);
		}

		/*
		 * update stress-like variable
		 */
		MutableDouble f = getVariable(e).getScalarVariable(VariableName.STRESS_LIKE_VARIABLE);
		f.setValue(computeSofteningLaw(wt.getValue()));

	}

	@Override
	public void FinalizeSolutionStep(Element e) {
		getVariable(e).save();
	}

	@Override
	public void CalculateOutputData(Element e) {
	}

}
